Probing Multi-Scale Energy Landscapes Using the String Method 
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A novel and powerful method is presented for the study of rare switching events in complex 
systems with multiscale energy landscapes. The method performs an umbrella sampling of the 
equilibrium distribution of the system in hyperplanes normal to the effective transition pathway, 
i.e. the transition pathway of a coarse-grained potential which need not to be known beforehand. 
Rather the effective transition pathway is determined on the fly in an adaptive way by evolving a 
smooth curve, called a string, that connects the initial and final regions and feels thermally averaged 
potential forces. Appropriate averages around the string determine the transition rates. Application 
to an example of solid-solid transformation of a condensed system is presented. 
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The evolution of complex systems often involves widely 
separated time scales. Well-known examples include nu- 
cleation events during phase transition, conformational 
changes of macromolecules, and chemical reactions. The 
separation of time scale is a consequence of the disparity 
between the effective thermal energy and typical energy 
barrier of the systems. The state of the system is confined 
for long periods of time in metastable regions in configu- 
ration space and only rarely switches from one region to 
another. The dynamics of the system effectively reduces 
to a Markov chain on the metastables regions. Find- 
ing the transition rates between the metastable states is 
a major computational challenge. For low dimensional 
systems, pathways and rates can be determined by iden- 
tifying the minimal energy paths (MEPs) which connect 
different minima on the energy surface by going through 
saddle points, and the curvature around these paths. For 
higher dimensional systems, such a procedure becomes 
impractical since the energy surfaces is typically dense 
in critical points and a multitude of dynamic paths con- 
tribute to the transition. The transition pathways have 
to be understood in an averaged sense as the superposi- 
tion of the many paths by which transition between two 
metastable regions in phase space takes place. 

Sophisticated numerical techniques have been devel- 
oped for finding transition pathways and rates |2|, 
though most methods apply when the energy landscape 
is relatively smooth (e.g. when transition is accomplished 
by a few isolated MEPs), and become less useful other- 
wise. A notable exception is the transition path sampling 
technique (TPS) || which applies to situations when the 
potentials are multiscaled. TPS samples the probability 
of dynamical trajectories conditional on their end points 
being in given regions of the configuration space. The 
transition rates are determined by reconstructing the un- 
conditional probability of the trajectories through um- 
brella sampling of trajectories whose end points are con- 
strained in overlapping sub-regions covering configura- 
tion space. TPS computes rates directly (i.e. without 
relating them to the geometry of the potential) . In this 



Letter, we propose an alternative approach which per- 
forms an umbrella sampling of the equilibrium distribu- 
tion of the system in hyperplanes normal to the transition 
pathways of a coarse-grained potential which need not be 
determined beforehand. Appropriate averages then de- 
termine the transition rates. 

We shall focus on the example of a system modeled by 



iq = -VV{q)+m 



(1) 



where 7 is the friction coefficient and £(t) is a white noise 
with (£j(i)£fc(0)) = 2jkBTSjkS(t). In molecular dynam- 
ics, (Q) is the high friction (or Smoluchowski) limit, but 
the method can be generalized easily to arbitrary fric- 
tion. We assume that the potential V(q) is multi-scaled 
in the sense that there is a gap between fc^T and the 
relevant energy barriers of the system. The system may 
certainly have other energy barriers that are comparable 
or smaller than fc^T. In this case, even though V may 
contain a very large number of local minima, the system 
experiences a much smoother energy landscape because 
of thermal effects and we are interested in the "thermally 
averaged" energy landscape: V := (V). The dynamics in 
both V and V reduce to the same Markov chain on suit- 
able metastable basins separated by the large energy bar- 
riers. In principle, one can study the transitions in such 
systems by first computing V and then applying meth- 
ods for smooth energy landscapes such as the nudged 
elastic band (NEB) method Q or the zero-temperature 
string method ||| to study the transition in the system 
with potential V. But in practice it is often impossi- 
ble to first find V. The motivation of the present paper 
is to develop a method that is based directly on V but 
computes the effective transition pathways and rates as- 
sociated with V. This is accomplished by a simple but 
important modification of the string method introduced 
in j3|] by incorporating thermal averages. The thermal 
average can either be the result of finite temperature, or 
the effect of noise introduced numerically to smooth out 
the small scale features of the energy landscape. 

It is useful to first review the zero-temperature string 
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method [Q. Let tp be a smooth curve, which we call 
a string, connecting two minima of V, A and -B. By 
definition, <y9 is a MEP if 



o=(v%)) 1 , 



(2) 



where (VV)" 1 " = VV — (VV ■ t)i and t is the unit tangent 
vector along tp, t = tp a /\tp a \. Equivalently tp is a curve 
which minimizes V in the hyperplane normal to itself. 
We will use suitable parametrizations of the string such 
that a — at A, a = 1 at B, One way of finding solutions 
of (H) is to follow the dynamics defined by: 



tp t = -(VV(tp)) ± +rl 



(3) 



where for convenience we renormalized time as t — > t/7. 
The scalar field r = r(a, t) is a Lagrange multiplier 
term added in (|^) to preserve some constraint on the 
parametrization of tp. For instance, a simple choice is to 
impose that tp be parametrized by normalized arclength. 
In this case (^) must be supplemented by the constraint 
(I Vol) a — which determines r. Other parametriza- 
tions - for instance by energy-weighted arclength - can 
be straightforwardly implemented by modifying this con- 
straint. See for details. 

This method is very effective if the energy landscape 
is smooth, as shown in W. In this case, it bears a lot 
of similarity with NEB But as we now show, it 

is the differences between the string method and NEB, 
namely the use of continuous curves with fixed and intrin- 
sic parametrization, that allows us to extend the string 
method to the case of rough energy surfaces. 

When the energy landscape is multiscaled, our ob- 
jective is to compute the effective transition pathways 
(MEPs in V) and rates, without computing first V. (Let 
us emphasis that since V is relatively smooth and the 
energy barriers associated with V are much larger than 
fcgT, we can assume that the effective MEPs are iso- 
lated.) To this end, we make an important modification 
of (J3j) , namely we replace the gradient of V at the right 
hand side of (j^) by an ensemble and/or temporal average 
of the actual force among a number of string configura- 
tions. We denote this ensemble by {tp^}, its mean by 
ip° := (ip u ), and therefore replace (E3) by 



-<vw)>- 



+ r°t° 



(4) 



where t° is the unit tangent vector along ip° and (■)~ L '° 
denotes the projection to the hyperplane normal to i°. 
To preserve parametrization of tp°, (^) must be supple- 
mented by some appropriate constraint which determines 
r°: for instance, the constraint (|Va|)a = corresponds 
to normalized arc-length. 

A practical way to create the ensemble {<^ w } is to in- 
troduce a finite temperature version of (||) through 



^=-(VF(^)) ± ' +r°t + «)- 



(5) 



where rf is a noise term added in order to simulate the 
effect of thermal averaging. We take rf to be a Gaussian 
process with covariance 

(^( ai t)<K,o))Jw (<) if " = "' (6) 

I otherwise, 

and compute all averages with respect to the statistics of 
jf . In particular, the mean of (||) is (|4|). 

The equilibrium density function for (^|) is given by 



t{q,a)=Z- 1 {a)e- 0V ^8s° {a) {q). 



(7) 



Here (3 = 1/fcgT, S°(a) is the hyperplane normal to 
tp°(a), 5s°( a )(o) is the Dirac distribution concentrated 
on S°(a), and 



Z(a) = / e-^d d q 

JS°(a) 



(8) 



is the partition function. (Q) can be determined by an- 
alyzing the Fokker Planck equation associated with (||). 
The mean string tp°(a) entering (Q) through S°(a) has 
to be determined self-consistently from 



tp°(a) 



q/i(q,a)d d q. 



(9) 



This equation implies that tp° is a MEP for V to leading 
order in fcsT. Indeed if V = V + SV with 8V comparable 
to ksT or smaller, the integral in (^|) can be evaluated 
by Laplace method and to leading order in fc^T ', only V 
contributes. This shows that tp° is the minimum of V(q) 
in the hyperplane perpendicular to tp° . 

(0) shows that the stochastic string performs an um- 
brella sampling of the equilibrium distribution of the sys- 
tem in the one-parameter family of hyperplanes normal 
to a MEP of V. (We stress again that this MEP needs 
not be known beforehand but rather is determined on the 
fly by the method.) Averages with respect to @) yield 
relevant quantities like, in particular, the transition rates 
between the basins of metastability visited by the mean 
string. This follows from a straightforward generalization 
of Kramers argument (see e.g. chap. 9 in Q) to multiscale 
energy landscapes. Let F(a) = — kgT In Z(a) be the free 
energy. Using the identity J d In Z/dada — lnZ(a), we 
obtain from (||) after integration by part 

F(a) = J ((t° ■ VV) ■ v °) a - t° a ■ tp))da. (10) 

Let cti, aii ■■■ be the successive values of a where F 
reaches a local minimum. Each otj corresponds to a basin 
of metastability visited by the mean string. Let ai2, «23j 
... be the successive values of a where F reaches a local 
maximum. Then 
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-/3F(a) 



tp° a \da 



tp° a \da, 
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a nontrivial test for reaction path algorithms [^) , and 5V 



as 



FIG. 1: Upper panel: the mean transition path tp° (black 
curve) on the disordered Mueller potential isolines. Also 
shown is the MEP, <p* , in V (red curve) obtained by inte- 
grating (^). Both curves coincide fairly well as expected. The 
region within the green curve is the root mean square dis- 
placement of the stochastic string indicating the typical size 
of the fluctuations of the string. Lower panel: the mean en- 
ergy along the stochastic string (black curve), (V((f)), com- 
pared with the smoothed energy along the MEP of V, V(tp*) 
(red curve). The blue curve is the total energy along <p° , 
V(<p°), which shows the large number of (irrelevant) minima 
and saddle points crossed during a transition. We took 21 
discretization points along the stochastic string and 20 real- 
izations. 



and so on. To leading order in fc^T, this expression is 
equivalent to 



\J F aa (ai)\F aa (a 12 )\ 
27r 7 K(ai)||^(ai 2 )f 



-PAFy. 



(11) 



where A_Fi2 := F{ol\2) — F(a±) is the free energy barrier 
from basin 1 to basin 2. 

In practice, (|B|) is solved by considering M realizations 
of the string, tp 3 , j = 1, . . . , M , and approximating the 
mean string as ip° = M^ 1 Yl^-i f 3 - After a number of 
steps, a reparamctrization step for (p° is applied. Once tp° 
has converged to its steady state value, no reparametriza- 
tion step is necessary anymore and, using ergodicity, one 
can supplement the ensemble average by a time average 
using <p° = (MT)- 1 Y^%L 1 fo<P j (t)dt to obtain better 
statistics. Other averages like ( |To[ ) are evaluated simi- 
larly. 

As a first illustrative example of the procedure, we 
consider a disordered version of Mueller potential. We 
take V to be the Mueller potential, originally invented as 
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5V(x,y)= ^2 (a k cos(2Tr(k 1 x + k 2 y)) 

klM = 5 +f3 k cos(27r(k 1 x-k 2 y))), 



[12) 



where a k ,Pk £ [— 1, 1]- With this choice |<5V| < 10 (com- 
pared to a main energy barrier of about 100 in V) and 
we shall assume that ksT = 5. The number of MEPs in 
V = V + SV joining the regions of metastability around 
the two deeper minima of V is of the order of 10 2 . This 
makes the complete sampling of the MEPs in V (say, 
using (||) with the full V, or NEB) impractical (and ir- 
relevant). The results of the stochastic string method 
are summarized in figure lj. The MEP in V is success- 
fully retrieved and, using ( 10J) , we also obtained the value 
AF12 = 109.0 for the free energy barrier between 1 and 
2, compared to the exact value of AF12 = 108.9 |t]]. 

Our next example illustrates the complex energy land- 
scapes for systems in condensed phases. We consider a 
crystal of 20 x 20 array of atoms interacting via a pairwise 
potential that has three minima at r% = l,r 2 = 1.1 and 
r3 = \/\ + l.l 2 . There are two basic equilibrium states 
for this systems, with the position of the (i, j)-th atom 
given by (r\i,r 2 j) and {r 2 i, r±j) respectively. Their rigid 
body rotations are considered to be equivalent states. 
The two states can be considered as the two variants of 
the martensitic phase. We are interested in the energy 
landscape that the system experiences when transform- 
ing from one variant to the other. 

For clarity we initiate the seed of transformation at 
the upper-left corner. Free boundary conditions are used 
in the simulation. The crystal transforms through the 
propagation of the twin boundary, oriented diagonally 
as shown in Figure 2. The twin boundary propagates 
through kink propagation along the twin boundary asso- 
ciated with the phase transformation of individual atoms. 
The energy landscape exhibits three scales: The largest 
scale associated with the position of the twin bound- 
ary, an intermediate scale associated with the propaga- 
tion of the twin boundary by one atomic distance, and 
a small scale associated with the kink propagation. For 
the stochastic string method, we choose fc^T to be in be- 
tween the energy barrier for kink propagation and twin 
boundary propagation. The energy landscape for the 
twin boundary to advance one atomic distance is shown 
in Figure 2. 40 points are used. In Figure 2 we also show 
the potential energy landscape associated with the MEP 
computed using the deterministic string method with 400 
points. Even though the details of the energy landscape 
is not resolved, as expected, the overall feature is cap- 
tured very well by the stochastic string method using 
only 40 points. 

In conclusion, we have introduced a powerful method 
to probe multiscale energy landscape and thereby deter- 
mine effective transition pathways, free energy barriers, 
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FIG. 2: Upper panel: Snapshot of the crystal during trans- 
formation. The color bar shows the scale of the local energy 
of each atom. Middle panel: The mean potential, (V(ip)}, ex- 
perienced by the crystal at finite temperature when the twin 
boundary moves by one atomic distance along the diagonal. 
This result was obtained using 40 discretization points along 
the strings. Also shown for comparison is the potential energy 
along one particular MEP (blue curve, at zero temperature) 
which shows the features of the potential associated with the 
transformation of individual atoms along the twin boundary. 
400 discretization points along the string were necessary to 
get this fully resolved result. Lower panel: The energy land- 
scape computed by the string method at finite temperature 
after the crystal is half transformed. Notice the appearance 
of three scales on the energy landscape. 



annealing schedule. This allows to further explore the en- 
ergy landscape and determine if more than one effective 
transition pathway is involved in a given transition. We 
are also investigating the possibility of using Metropo- 
lis Monte-Carlo scheme instead of the Langevin equa- 
tion (j|) to compute averages with respect to the equi- 
librium distribution in (Q). This speeds up convergence 
since it requires to compute V only instead of its gradi- 
ent. Finally, in its present form the method performs the 
coarse-graining in an implicit way, i.e. no reaction coor- 
dinate has to be known beforehand. In certain systems 
where reaction coordinates can be easily determined, the 
method can be applied to these coarse-grained variables 
alone, which might dramatically speed up convergence if 
the reaction coordinates are low dimensional. 

We thank David Chandler and Bob Kohn for helpful 
discussions. The work of E is supported in part by NSF 
grant DMS01-30107. 
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and transition rates in complex systems. The method 
performs umbrella sampling in the hyperplanes perpen- 
dicular to an effective transition pathway which is deter- 
mined on the fly in an adaptive way. The method can be 
further improved in various ways. For instance the tem- 
perature of the string can be periodically changed on an 



